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Recent developments in the search for topological superconductivity have brought lattices of 
magnetic adatoms on a superconductor into intense focus. In this work we will study ferromagnetic 
chains of adatoms on superconducting surfaces with Rashba spin-orbit coupling. Generalising the 
deep-impurity approach employed extensively in previous works to arbitrary subgap energies, we 
formulate the theory of the subgap spectrum as a nonlinear matrix eigenvalue problem. We obtain 
an essentially analytical description of the subgap spectrum, allowing an efficient study of the 
topological properties. Employing a flat-band Hamiltonian sharing the topological properties of 
the chain, we evaluate the Z-valued winding number and discover five distinct topological phases. 

Our results also confirm that the topological band formation does not require the decoupled Shiba 
energies to be fine-tuned to the gap centre. We also study the properties of Majorana bound states 
in the system. 

PACS numbers: 73.63.Nm,74.50.-l-r,74.78.Na,74.78.Fk 


I. INTRODUCTION 

The physics of magnetic impurities in superconduc¬ 
tors has been studied extensively since the seminal works 
by Yu, Shiba and Rusinov.EH^ A magnetic moment lo¬ 
cally disrupts the superconducting condensate, giving 
rise to subgap bound states localised in the vicinity 
of the magnetic moment. In the past two decades. 
Scanning Tunneling Microscopy (STM) studies of mag¬ 
netic impurity ph ysics have complemented this theoret¬ 
ical understanding.®^ The interest in these systems was 
renewed after proposal^®i^to realise ID and 2D topolog¬ 
ical superconductors by lattices of magnetic adatoms on 
superconducting surfaces. The promising observatior^ 
of signatures of Majorana bound states (MBSs) in fer¬ 
romagnetic chains has brought the topic in the intense 
focus recently. Similar signatures are also claimed to 
have been observed in another recent experiment .^Mag¬ 
netic chains, with their respective advantages and weak¬ 
nesses, offer an interesting alternative to the nanowire 
systemJiSlti^ as a route towards topological superconduc¬ 
tivity. 

In this work we concentrate on the properties of ID fer¬ 
romagnetic adatom lattices deposited on a 2D supercon¬ 
ducting film with a Rashba spin-orbit coupling (SOC). 
We will study the situation where each magnetic moment 
binds a single subgap Shiba state with energy A , 
where the dimensionless coupling a = JSvqtt is deter¬ 
mined by microscopic parameters - vq is the density of 
states, J is the exchange coupling and S is the magni¬ 
tude of the impurity spin - and A is the pairing gap 
of the superconductor.®i21 xhe Shiba states are concen¬ 
trated around the magnetic moments with wavefunctions 

-r/€ -r/€ 

decaying as ^ 372 “ in 2D and — in 3D superconduc¬ 
tors, where the exponential decay is controlled by the su¬ 
perconducting coherence length ^o- The wavefunctions 
are more slowly decaying in 2D, which was qualitatively 



Figure 1. Schematic figure of the ferromagnetic chain of atoms 
embedded on a two-dimensional superconductor. 

confirmed in a recent experiment.^ The slow spatial de¬ 
cay of the bound states at length scales below ^ has 
important physical implications. The separation of the 
magnetic moments can easily be much shorter than the 
coherence length, thus allowing a significant hybridisa¬ 
tion of bound states separated by dozens of neighbours. 
Therefore effective theories of subgap bands generically 
exhibit a long-range hopping, adding its special proper¬ 
ties to the problem. The long-range hopping nature of 
the Shiba states will give rise to rich topological proper¬ 
ties, especially in 2D magnetic lattices where it is possible 
to achieve dozens of dif ferent phases and Chern numbers 
much larger than unity.l^^*^ 

The topologi cal properties of magne tic chains with 
helicaP® 2 lllH^ and ferromagneticPSESI order have al¬ 
ready been studied in considerable detail. An interesting 
recent step towards a more realistic description of Shiba 
systems was taken by the generalisation of the theory of 
topological Shiba bands to the multichannel case to ac¬ 
commodate multiple bound states on a single impurity.l^ 
That work, like most of the studies capturing the micro¬ 
scopic structure of the Shiba states, is restricted to the 
deep-dilute impurity limit where the energy of the de¬ 
coupled Shiba states reside close to the gap centre, i.e. 
a « 1. This is motivated partly by the signihcant theo¬ 
retical simplification arising from the linearisation of the 
spectral problem in energy, which allows an elegant and 
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simple formulation in terms of an effective HamiltonianP 
Our work takes a step towards a more general descrip¬ 
tion by formulating a theory of Shiba bands beyond the 
deep-impurity regime. This is important, since tuning 
the parameters controlling the value of a is difficult in 
experimental setups. As highlighted in our work, the mi¬ 
croscopic parameters do not need to be fine tuned to the 
deep-impurity regime to reach the topological phase. 

In Sec.|nJ we formulate the subgap spectral problem of 
a ferromagnetic Shiba chain as a nonlinear matrix eigen¬ 
value problem and obtain an essentially analytical solu¬ 
tion for the dispersion and wavefunctions valid for ar¬ 
bitrary subgap energies. In Sec. m we show that the 
ID chain has an emergent chiral symmetry which en¬ 
ables us to defining a Z-valued winding number invariant 
M classifying the different phases. In Sec. |IV[ we dis¬ 
cover five distinct topological phases with winding num¬ 
bers M = —2, —1,0,1, 2, discuss the implications of the 
observation and compare the results to the deep-impurity 
formulation. We study the MBSs supported by the chain 
in Sec. |V) extracting the spatial decay of their wavefunc¬ 
tions, topological gaps and energy splitting. In Sec. |VI[ 
we summarise our results and discuss their implications. 


II. MODEL AND THE SPECTRAL PROBLEM 


A. Subgap spectrum as a nonlinear eigenvalue 
problem 


The starting point of our study is the system de¬ 
picted in Fig. consisting of a 2D superconductor sur¬ 
face decorated by magnetic impurities arranged in a chain 
with lattice constant a. The system is described by the 
Bogoliubov-de Gennes Hamiltonian 


n = 



M + otji{kyax 



+Atx - J ^(S • (Ti)S{r - Ti). 

i 


( 1 ) 


The matrix structure of Eq. Q corresponds to the 
Nambu basis 'k = ('i/'ti''“i ai cor¬ 
respond to Pauli spin matrices in the particle-hole and 
spin subspaces, respectively. Here k'^/2m — /r is the ki¬ 
netic energy of the electrons, an is the spin-orbit coupling 
and A describes the pairing amplitude of Cooper pairs. 
The vector r is the position of the electron, whereas 
describes the impurity positions. We will focus on the 
case where all magnetic moments point in the z direc tion. 
Consequently the system is in the symmetry class 
Inserting the Hamiltonian density into the Bogoliubov-de 
Gennes equation = £’![', we obtain 

[E - (^k + aR{kyax - k^Oy)) - At^] ^'(r) 


where we have introduced = fc^/(2m) — /r. We can 
make further progress by restricting the system to a one¬ 
dimensional chain in the x direction, yielding, as seen in 
Appendix A, the equation 


E + Tx 

I — a , - Sa 


VA2 - A2 






( 3 ) 


where a = JSvqtt = ^JSm, and Jr is essentially given 
by the propagator of the 2D bulk electrons. The specific 
form of Jb in this case was first evaluated by Brydon et 
and is presented in Appendix A. In the deep-dilute 
limit this equation can be linearised in E and \/y/kpa, 
where kp is the Fermi wave number and a is the chain lat¬ 
tice constant; this allows projection onto the low-energy 
single impuriW bands, resulting in an effective two-band 
Hamiltonian.lSl However, in this work we are interested in 
the behaviour of the system for arbitrary subgap energies, 
thus requiring a more general approach to the problem. 
As outlined in Appendix A, it is convenient to work in 
the basis of the eigenstates of Txa^, in which case we 
can obtain a nonlinear eigenvalue problem (NLEVP) for 
A = (A -b E)/yJA‘^ - and 


(aX^-^X 

bX cX? 

-bX 

4A — a —Xd 

—cX^ 

—Xd aX^ + xA 

\ -Xd 

—c bX 


-Ad \ 
c 

-bX 

--A — a) 

a / 


( 4 ) 


Here a, b, c, d are N x N matrices (where N is the length 
of the chain), and similarly A is shorthand for AIjvxAf- 
The coefficient matrices are of the form 


bij — 


\/A2 — E^ 

- ^ -(4" 


dj — z 


.V'A2- 


( 5 ) 


2m 


{I4 i^ij) I4 


dij — 


2^ (A i^ij) d- l/' {Xij)), 


where = Xi — Xj denotes the difference between po¬ 
sitions of magnetic moments and are given by the 
special function expressions presented in Appendix A. It 
is important to note that the submatrices have a nontriv¬ 
ial energy dependence through the energy dependent co¬ 
herence length = vp/^/A'^ — E"^ (not to be confused 
with the kinetic energy ^k), which in principle compli¬ 
cates the solution considerably. In the limit where the 
energy-independent coherence length = vp/A goes to 
infinity this nontrivial dependence vanishes. The prob¬ 
lem Q is a representative of NLEVPs A{E)'^ — 0, where 
A(E) is a matrix-valued non-linear function of E and 
can be thought of as a generalisation of the usual linear 
Schrodinger problem where A{E) = H — El. 
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On physical grounds we know that since we are deal¬ 
ing with a system of N magnetic moments with a sin¬ 
gle s-channel bound state, Eq. has N positive and N 
negative energy solutions describing the magnetic subgap 
band. Most of the previous works solve the problem in 
the deep-impurity limit by linearising the model, yielding 
a linear problem = £'4' with an effective 2N x 2N 

Hamiltonian 

= ( 6 ) 

\^—^ac' aa — lij ' ' 

where d = lim£;_>oa, c = lim£)_>o c. This description is 
valid for deep impurities awl close to the gap centre 
coupled by a weak hopping 1/^/kpa ^ 1. One of our 
central results presented below is that by relaxing the 
conditions leading to the two-band model we can show 
that the topologically nontrivial phases extend outside 
the deep-impurity regime. Therefore it is not necessary 
for the decoupled Shiba states to lie close to the gap cen¬ 
tre to achieve topological phases, and there is no need to 


tune their energies if the hybridisation, which can be con¬ 
trolled by the distance between the magnetic moments, 
is sufficiently strong. 


B. Solution to the nonlinear eigenvalue problem 

In the case of a periodic or infinite system, the problem 
can be vastly simplihed by going over to reciprocal space. 
The submatrices a, 6 , c, d are all translationally invariant, 
allowing a Fourier transform 

= (7) 

which reduces the equation to a 4 x 4 NLEVP. In this 
form we can hnd an expression for the bulk spectrum 
from the usual eigenvalue condition by requiring that the 
determinant of the matrix in Eq. Q vanishes. This gives 
us a transcendental equation for the energy. 


E(k) = ±A. 


'_ (gfc + ~ ~ 4(afc^fc + Ckdk)'^ _ 

(®fc + + Cfc + - l/a)2 - 4(afe6fc -|- CkdkY + 4:{al + cl) ' 


In comparison, the expression for the energy of the lin- are, up to a normalisation factor, 
earised two-band model is 


( 8 ) 


E‘^^{k) = ±A^J (1 - adk)'^ + a'^cl (9) 

- note that this is not a transcendental equation, be¬ 
cause E is set to 0 in the expressions for at and Ck and 
hence it gives E directly for any coherence length. How¬ 
ever, Eq. ([^ can be solved numerically for any coherence 
length with standard methods for transcendental equa¬ 
tions. As shown in Appendix B, it turns out that replac¬ 
ing with ^0 in the four-band model gives an excellent 
approximation of the energy. This has two important 
consequences. Firstly, after setting ^e —t Co on the right- 
hand side, Eq. Q essentially provides an explicit solution 
of the subgap energy bands typically within accuracy of 
the order of 10“^A or even much better. The accuracy 
of the solution could be numerically improved systemat¬ 
ically by iteration as also shown in Appendix B, but for 
all practical purposes the zeroth order solution is suffi¬ 
cient. Secondly, since the replacement Cb Co in the 
bulk spectral problem is seen to lead to negligible correc¬ 
tions to the energies, we also employ this approximation 
when we solve the NLEVP in real space to study the 
MBSs. 

Remarkably, the eigenspinor can also be solved ana¬ 
lytically from the 4x4 NLEVP in reciprocal space. The 
components of the eigenspinor '^\{k) = (xi, X 2 , X 3 , X 4 )^ 


Xi = —A 


c" - d" + 


2 X^d^ — {ad — 6 c)^'' 


X 3 = A(ac — bd) — c-\- 


A^ — (a^ -I- c^) 

X{ab -I- cd -I- Xb){bc — ad + Ad) 
A^ — (a^ -I- c^) 


A(a6-I-cd-I-A6) X{ad—bc+\d) 

^2 ^ r, } 7) I Q., “t“ E7) Q I Qx ^3 


A^ — (a^ -I- c^) 


A^ — (a^ + c^) 


X 4 = - [Xbxi -I- (a — A)x 2 -I- Adx 3 ]. 
c 


( 10 ) 


where we have suppressed the k index on the functions 
a, b, c, d. The eigenvector |£+) corresponding to the pos¬ 
itive energy can be found by inserting the appropriate 
solution E{k) > 0 of Eq. (|^ into the expressions (10) 
for the components. Importantly, the negative energy 
solution \E-) corresponding to E{k) < 0 is related to 
the positive energy solution (at the same k) through 
C|£+) = |£-) where C = TyCjy in the original basis of 
Eq. ([^ , ov Gz® Ux in the rotated basis in which solution 
( [I^ is obtained. The expression Gz ® stands for a 
4x4 matrix where the Pauli matrices do not refer to the 
electron spin anymore. 
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III. TOPOLOGICAL CONSIDERATIONS 


A. ID flat-band Hamiltonian with chiral symmetry 


In this section we will present the topological prop¬ 
erties of the subgap Shiba bands described by Eqs. Q 
and (10). In the periodic table of topological insulators 
and superconductors, d-dimensional gapped systems are 
typically classified by studying associated d-dimensional 
(Bloch) Hamiltonians and their symmetries. However, 
while we started with a 2D Hamiltonian ([^ we obtained 
subgap energy bands that are essentially localised in the 
vicinity of the ID chain of magnetic atoms. In addition, 
the eigenstates \E±) are parametrised by the ID wave 
vector k and so the subgap spectrum is effectively one¬ 
dimensional. In the deep-dilute limit it is possible to 
reduce the NLEVP to a linear problem in terms of an 
effective Hamiltonian of type ® and study the topolog¬ 
ical properties by treating H^"as a hona fide ID Hamil¬ 
tonian. However, in our case another strategy must be 
adopted due to the nonlinear formulation of the problem. 
Therefore we define an effective Hamiltonian in terms of 
the projectors of the subgap bands as 


H=Y^ E,\E,){E,\ = ^ E,P, = E+P+ + E_P_. 

i/=± i/=± 

( 11 ) 

The low-energy properties of the initial 2D system coin¬ 
cide with those of the ID Hamiltonian 0- When we are 
only interested in the topological properties of the sys¬ 
tem, the precise details of the energy bands are irrelevant 
and with no loss of generality we can adiabatically flatten 
them E±{k) —>■ ±I and study P = P+ — P-. Notice that 
the property C\Ef) = \E-) implies the anticommutation 
relation {C,H} = 0, i.e., the ID system has the chiral 
symmetry and belongs to the symmetry class BDplJ^ 
as the original Hamiltonian ([^ if motion is restricted to 
the X direction {ky = 0). 


B. Winding number 


In noninteracting ID systems with chiral symmetry 
the topological classification is Z-valued, thus the system 
can support an arbitrary number of different topological 
states. These states can be most conveniently obtained 
by evaluating the winding number invariant given in the 
form 

1 r'^/a- 

Af=— dktT \CH-^dkH] . (12) 

diri 


where C = = C ^ is the chiral symmetry operator 

and iJ is a fully gapped ID Bloch Hamiltonian! ^^ ! ^^ ^ As 
discussed above, in the transformed basis the chiral sym¬ 
metry operator for this system is C = (72 0 o'-.-. After 
inserting the flat-band Hamiltonian H in Eq. (12) and a 
few lines of algebra, outlined in Appendix C, the winding 


number can be expressed in terms of the band projectors 
as 




1 

iri 



dk tr 


CP+dkP+ 


(13) 


In the next section, we will analyse the topological phase 
diagram by evaluating (13) as a function of the system 
parameters. 


C. Z2 invariant 


The existence of chiral symmetry is typically ap¬ 
proximate in real systems and topological classification 
based on particle-hole redundancy that is built into the 
Bogoliubov-de Gennes formalism is expected to be less 
sensitive to disorder. Therefore we also consider Kitaev’s 
Z 2 invariant measuring the ground state fermion parity. 
The phase transition between different parities is associ¬ 
ated by the energy gap closing E{k) = 0 a,t k = 0,7r/a. 
This is useful, as the phase boundary between of differ¬ 
ent Z 2 phases can be explicitly obtained from the NLEVP 
matrix in reciprocal space, as discussed in Ref. [331 We 
are especially interested in the phase diagram as a func¬ 
tion of the single-impurity Shiba energy determined by a 
and the hybridisation parameter kpci. Setting E{k) = 0 
at k = 0,7r/a, we can solve the gap closing lines as a 
function of a and kpa as 



(14) 


This equation defines the phase boundary between re¬ 
gions of different parity of the Z 2 invariant, and is valid 
for arbitrary coherence lengths. The corresponding equa¬ 
tion for the two-band model is given by 


1 


a = — 

Ofc 


fe=0,7r/a 


(15) 


allowing for simple comparison of the Z 2 phase diagram 
of the exact four-band solution and the two-band approx¬ 
imation valid in the deep-impurity limit. In the presence 
of chiral symmetry the Z 2 phases are captured by the 
winding number parity so considering it separately brings 
little new information. However, formulas (14) and (15) 
yield semianalytic phase boundaries of the different par¬ 
ity phases and make their analysis very convenient with¬ 
out the need to perform integral (13) everywhere. 


IV. TOPOLOGICAL PHASE DIAGRAMS 


In the previous section we derived analytical formulas 
describing the FM Shiba chain. In this section we will 
use those results to examine the topological properties of 


the chain. To that end we have used Eqs. (13) and (14) 


to calculate the topological phase diagram of t^ 


le system. 
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Figure 2. a) Topological diagram of the four-band model, 
created by plotting the winding number as a function of kpa 
and a. The parameters used are = 50a, ? = 0.1. Winding 
numbers from —2 to +2 are represented. As kpa increases 
the winding number configurations (quasi-)periodically flip 
signs; this happens more often with increasing <;. b) Same 
parameters as in a), but focused on a narrower range of kpa 
values, c) Energy of the lowest-lying positive state for an 
infinite system as a function of kpa and a. Parameters used 
are the same as in a). The topological phase transitions are 
seen here as bulk gap closings. 


As seen in Fig. the phase diagram thus obtained 
reveals the presence of several topological phases, in ac¬ 
cordance with the BDI classification of the FM chain. 
In total we found five phases, corresponding to wind¬ 
ing numbers 0, ±1,±2. This feature is present in both 
the general model and in the linearised two-band model, 
a feature missed in previous works which either consid¬ 
ered the parity of the winding numbeili^ or examined 
a parameter area too narrow to contain all the wind¬ 
ing numbers.^ Often, nontrivial regions appear for pa¬ 
rameter values where the validity of the two-band model 
breaks down, so a reliable identification of the phase dia¬ 
gram requires the full machinery of the four-band model 
derived in the present work. As seen in Fig. [^c), the 


topological phase transitions are accompanied by closings 
of the gap as expected. A phase with winding number 
N supports lA/”! degenerate MBSs; further, domain walls 
between regions with different winding numbers Mi and 
M 2 can support lA/*! — A 2 I bound states. However, while 
the energy gaps for phases lA/”! = 1 can be as high as 
0.25A, the jA/"! = 2 phases generally have smaller gaps; 
this naturally provides a complication in realising domain 
walls with more than two MBSs. The general behaviour 
of the energy gap will be discussed in more detail in the 
next section. Notice that significant regions of the non¬ 
trivial topological phases extend to values a > 1.5 and 
a < 0.6 for small values of kpa (large hybridisation), 
translating to bound state energies between —0.4A and 
0.5A. Thus, provided that the separation of the moments 
can be tuned in the fabrication of the chain, there is no 
need for tuning the bound state energies close to the gap 
centre. 


It is interesting to compare the accuracy of the two- 
band mode l to the exact four-band solution. Several re¬ 
cent paper^SMOni have considered a similar chain with 
a helical spin configuration in an intrinsically 3D sys¬ 
tem rather than the ferromagnetic planar version with 
Rashba SOC that is the focus of this article. As seen in 
Ref. 1331 the two-band and four-band models are in fairly 
good agreement in the helical system; however, in the fer¬ 
romagnetic case the differences of the models are more 
pronounced. While the two-band and four-band models 
indeed agree in the deep-dilute limit a « 1, kpa ^ 1, 
the convergence is much slower than in the helical 3D 
model. In order to compare the two models, we plot the 
Z 2 phases of the four- and two-band models, using equa¬ 
tions (141 and in Fig. 1^ Even though the two-band 
model captures the behaviour reasonably in the dilute 
regime kpa > IOtt, it is evident that non-negligible de¬ 
partures remain. In the dense regime kpa ^ 47r the dif¬ 
ferences become qualitatively significant and the deep- 
impurity approximation starts to break down. 


Having analysed the impact of the parameter a, we 
move on to consider the behaviour of the topological 
phases as a function of the normalised SOC strength 
g = mapi/kp. In Fig. [^we plot the winding number as a 
function of kpa and for different values of a. While the 
trivial regions do grow when the single-impurity bound 
state is moved away from the gap centre a = 1, for the 
parameters chosen here, the M = ±2 phases grow as well; 
hence, lowering a in this case makes even parity more 
likely. While the linearised model may not quantitatively 
agree with the four-band model, its phase diagram shows 
many of the same characteristics. To make contact with 
previous work, we present in Fig.|^c) a topological phase 
diagram similar to that of Fig. 1 in Ref. [T^l The evalua¬ 
tion of the winding number reveals a richer structure of 
topological regions as well as entirely new non-trivial re¬ 
gions compared to the Z 2 classification presented in that 
work. 
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Figure 3. Z 2 phase diagram of the exact four-band model 
(red) and the two-band deep-impurity (blue) model. Param¬ 
eters used are ^0 = 50a, ? = 0.01. Note that the agreement 
is not accurate even in the deep-impurity a = 1 regime. The 
agreement improves as kpa increases and the system becomes 
more dilute. 


V. MAJORANA BOUND STATES 


In the previous section we analysed the topological 
properties of the system, finding values for the winding 
number between —2 and -1-2. In practice, realisation of 
these phases may be limited by factors such as system 
size and temperature. In this section we will study the 
energy scales and spatial decay of the Majorana bound 
states in finite chains with open boundary conditions. 
This is carried out by numerically solving the 4A^ x 4A^ 
NLEVP (|^ in real space. 

A robust topological phase requires the existence of a 
well-defined energy gap between the MBSs at zero energy 
and the bulk states. In the limit of an infinite system, this 
gap can be straightforwardly calculated using Eq. In 
Fig. a) we have plotted the energy gap of the system 
as a function oi kpa and a. Excepting an area of low 
gap size around certain k = 0 transitions, the energy gap 
is generally of the order of O.IA, reaching values around 
0.25A at the center of the \Af \ = 1 region; assuming a gap 
size similar to the 1.36 meV seen in the experiment on 
Pb surface in Ref. this corresponds to temperatures 
« 4 K. While the parameters for a realistic system may 
not reach this maximum it is clear that there are wide 
regions where the gap size is large compared to temper¬ 
atures that can realistically be reached in experimental 
setups. The \Af \ = 2 phases generally have a much lower 
bulk gap, often by an order of magnitude. For finite 
systems the energy gap and the MBS splitting will de¬ 
pend on the chain length. In Fig. b) and c) we have 
plotted the gap and MBS energy as a function of chain 
length in the lA/”! = 1 and lA/”! = 2 region, respectively. 
In the lA/”! = 1 region the dependence of the energy gap 
on the chain length is negligible, and for relatively short 
chains with N < 50 the MBS spitting is two orders of 




Figure 4. Topological phase diagram plotting the winding 
number as a function of kpa and ?. The color scheme is 
the same as Fig. a) Parameters used are a = 1, ^0 = 
50a b) Same diagram, but a = 0.85. The trivial regions 
at the center and edges of the figure are larger than in a), 
but so are the \J\f\ = 2 phases, even at the expense of some 
trivial regions, c) Topological phase diagram of the two-band 
model. Parameters used are a = 1, ^0 = 5a; compare Fig. 1 
in Ref. [121 which plots the parity of the invariant for the same 
parameters. 


magnitude smaller than the energy gap. The oscillating 
splitting of the MBSs localised at each end of the wire 
will first go down exponentially after which it settles on 
a slower algebraic suppression as discussed previously for 
wavefunctions in Ref. [^ The algebraic tail of the MBS 
splitting could, in principle, be harmful for quantum in¬ 
formation applications, though the absolute value of the 
splitting is orders of magnitude smaller than the gap. For 
the lA/”! = 2 phase the energy gap is first reduced by an 
increasing chain length before converging to a constant 
value. Observation of the double MBSs at each end is 
limited by the large splitting seen in short chains, hence 
requiring larger system sizes and lower temperatures. 

Besides the splitting behaviour we are also interested 
in the spatial decay of the MBS wavefunctions. This 
was previously analysed in the case of a helical chain by 
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Figure 5. a) Bulk energy gap as a function of kpa and a. Parameters used are ? = 0.01 and = 50a. The black lines 
correspond to topological phase transitions; the wide black zones are areas of low bulk gap. The corresponding winding number 
diagram is illustrated in Fig. |^in Appendix B. b) Dependence of bulk and MBS energy on the length of the chain N in the 
Af — 1 phase. Parameters used are kpa = 20, a = 1, c = 0.01, = 50a. c) Dependence of bulk and MBS energy on the length 

of the chain in the N = 2 phase. Parameters used are kpa = S.Ttt, a = 0.8, ? = 0.01, = 50a. 



Figure 6. Spatial decay of Majorana wavefunctions. a) Wave- 
function in A/” = 1 phase. Parameters used are N = 500, 
kpa = 20, a = 1, ? = 0.01, = oo. The fitted function 

is /(n) « 0.43e“°'^®" -f 6.84n"Mn-2(n/0.0011). b) Wave- 
function in M = 2 phase. Parameters used are N = 800, 
kpa = 17.74, a = 0.73, ? = 0.01, = oo. The fitted function 

is /(n) « 0.31e“° °^®” -f 1.01n-° '^® ln-2(n/0.07). 

Pientka et An analytical solution of this problem is 
beyond the scope of this work, but numerical parameter 
fitting indicates that the envelope of the wavefunctions 
decays in a similar manner as those found in the helical 
chain. As seen in Fig. the decay is exponential over 
short distances but over longer distances (but still < ^o) 
takes a character of a;“^ ln(a;/xo)“^, where /3, 7 S K-i-. 
This decay of the wavefunctions reflect the behaviour of 
the MBS energy splitting which has the same origin. In 
general, the MBSs seem to be more localised in the lA/”! = 
1 phases, while the lA/”! = 2 phase features stronger finite- 
size effects which are also apparent in Fig. 

VI. DISCUSSION AND OUTLOOK 

Motivated by recent experimental developments in 
topological chains, we studied the case of a ferromag¬ 
netic chain of adatoms embedded on a superconductor 
with Rashba spin-orbit coupling. Building upon previ¬ 
ous works limited to the deep magnetic impurity regime, 
we formulated a theory that can address ID Shiba bands 
at arbitrary subgap energies. This allowed us to make 
a number of novel observations regarding the topological 


states in the chain. 

We found that the system, in accordance with its BDI 
classification which allows a Z-valued topological invari¬ 
ant, supports five different topological phases, some of 
which have hitherto not been reported in this setup. 
These phases may support zero, one or two Majorana 
end states and their interfaces may support up to four 
Majorana bound states. Our analysis revealed that en¬ 
ergy gaps of the most robust phases may be as high as 
0.2A and that the decoupled impurity energy does not 
need to be fine-tuned close to the gap centre for obt aining 
robust topological phases. In the dense-chain limilpSlMl^ 
where the orbitals of the magnetic atoms have a direct 
overlap, the mechanisms of topological phases are quali¬ 
tatively quite different, bearing a strong resemblance to 
the proximity coupled nanowire physics. However, it is 
remarkable that robust topological gaps can be obtained 
also in the studied Shiba limit with bound state energies 
well away from the gap centre. The sources and effects of 
disorder in the long-range hopping model are expected to 
be quite different from the conventional nanowire setting 
and are left for future studies. Considering the richness 
of the topological properties, and the fact that the system 
allows for local probing in STM, the studied system is one 
of the most promising platforms for studying topological 
superconductivity. 

The main ob stacle for applications of topological 
superconductivitjlSSMI may be the slow spatial decay of 
the edge modes, which results in an algebraic energy 
splitting of the Majorana end states. However, already 
in relatively short wires the splitting could be orders of 
magnitude smaller than the bulk energy gap, as found in 
our numerical calculations. Our results are generally en¬ 
couraging for the prospect of an experimental observation 
of topological superconductivity in the studied setup. 
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Appendix A: Derivation of the NLEVP 


Multiplying the equation by the inverse of the matrix on 
the left-hand side, and transforming back to real space 
In this appendix we will derive the nonlinear eigenvalue yields 
problem for the four-band Rashba model. Our starting 
point is Eq. (§ 

[E - (Ck + anikyaj; - ka^dy)) - At^] ^'(r) 

=-J^S'cr^(5(r-rj)4'(rj> (^-1) 

3 


which we then Fourier transform, giving us 

[E - (^k + a.R{ky(Jx - kx<7y)) Tz - Atj,] 

= (A-2) 

j _ I 


(1 -k 


with the identification 


F! JE{xi - Xj)Sa^'S{xj). 


i¥^i 

(A.3) 


/ dV 

- (^k + anikyd^ - k^dy))T^ - At^\' 


(A.4) 


By defining 


calculate, namely 


M± = 


E -\- ^±Tz -\- ATx 
A 



(A.5) 


where k = yk'^ + k'^ and = ^k ± ctuk, we split Je 
into two helicities corresponding to each sign, 


/ dk 

(A.6) 


Itir) 

Itir) 

liir) 

liir) 


N± f 
(27r)2 io 
N± f 
(27r)2 Jo 
N± f 
(27r)2 Jo 
N± f 
(27r)2 Jo 



dC 

d^ 

d^ 


cos 0 

E^-S3- A2 


^^ie+ik±{i) 


r cos 6 


A2 - ^2 _ ^2 


(^)r cos 9 

E^-S3- A2 


^i9-\-ik± (O’" cos 9 


A2 - ^2 _ ^2 ■ 


(A.7) 


Here we have introduced N± = 1 =p <;/\/\ + <^2 and 
k±[i) = /cf(vVV? T <?) + C/'I'fVTV?, with the nor¬ 
malised spin-orbit coupling <r = mau/kp. To make fur¬ 
ther progress, we restrict our attention to subgap ener¬ 
gies, i.e. E < A. For r = 0, all integrals but the third 
vanish trivially, and is simplified enough to be cal¬ 
culated using standard contour integral methods. For 
r > 0 we must employ Bessel and Struve functions. Set¬ 
ting the chain to be along the x-direction, we can now 
rewrite Eq. ( |A.3[ ) into the form of Eq. (|^ with 


We assume that the wavevector k is centered around the 
point « 0, or, equivalently, in the neighbourhood of 
the Fermi surface for both helicities. In spherical coor¬ 
dinates, we can then substitute k for in the integrals 
as appropriate. This gives us four different integrals to 


Je{x) = f [(/i (x) -k It {x))tz + (/a (x) - It {.x))Tzdy\ 

-kf (E -k Atx) [(/g” (x) -k It {x)) + (/^ (x) - /4^(x))cry], 

(A.8) 

where 
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lf(r) = Af±Im [Jo((fcF.± + i^E^)\x\)) + iHQ{{kF,± + *?£;^)|a;|))] 

iti’’’) = -iA^±sgn(2;)Re [iJi((/cF,± +*?E^)|a^l)) + H_i{{kF,± +iC£;^)kl))] 

[U{kF.,± + *eF')kl)) + ^M{kF,± + *eF')|2;|))] 

ktir) = [iJi((A:F,± + *^F^)|a;|)) + i?-i((fcF.± + iCF^)|3^l))] • 


Here we have used the previously mentioned Bessel (Jo,i) 
and Struve (Hq.-i) functions along with the two Fermi 
momenta = kF^s/l + T ‘t) and the quantity 

= Co/\/l — E'^jIS?. The problem can be further 
simplified by working in the basis of the eigenstates of 
'brO’zj which correspond to the eigenstates of the single¬ 
impurity problem. With the shorthand \tx<7z)i where for 
example |-|- f) = |+)r^ ® | the spinor in the new 

basis is 


'f'. = ((+ti^,) (-il'i',) i+mj) (A.io) 


In the transformed basis, we obtain a AN x 4iV nonlinear 
eigenvalue problem for a chain of length N. By defining 
the matrices 


_ J^2 

Qij = {I^ i^ij) + I3 + ^ij 

bij = i^ij) ~ ^2 


Ci -i — t 


VA 2 -f ;2 _ 


(A.11) 


2 m 


(^4 {,'^ij') li 


d'ij — 2 rrS^^ 


we can 


transform Eq. (|^ into the NLEVP in Eq. Q. 


Appendix B: Validity of the energy-independent 
coherence length approximation 

The purpose of this appendix is to give some motiva¬ 
tion as to why using instead of ^f in the four-band 
model is an excellent approximation. We first note that 
for a given we have A Coi so it is evident that the 
real energy must lie within the interval [/inf, /sup], where 
/inf /sup stands for the infimum/supremum in the interval 
\ G [Co,oo) for the function / defined by E = A/(^f), 
i.e. the right-hand side of the expression for the energy in 
Eq. ([^. This, however, does not narrow the error down 
too much as f[C) can vary greatly between ^ = /o and 
^ = oo. 

One way to approach the problem is through the use 
of contractions. A contraction 5 : K —)■ K is a Lipschitz 
continuous function, whose Lipschitz constant Kg < 1. 
In other words, if g is a contraction and K is equipped 
with the usual metric, we have 


(B.l) 




-^-£ 2 ] 


10 '’® 


0 


0.1 


0.2 0.3 0.4 

l/?o 


Figure 7. a) Energy as a function of k for the first approxima¬ 
tion £0 as well as the numerically exact fixed-point solution e 
of Eq. I^. Parameter values are kpa = 20, 0 = 1,? = 0.01, 
^0 = 50a. Iterations are omitted here as they would be in¬ 
distinguishable from the fixed-point solution, b) Difference in 
energy between the solution of the transcendental equation 
and several iterations of the solution for £„, plotted as a func¬ 
tion of the inverse coherence length at zero energy. Parame¬ 
ters used are kpa — 20, a = 1, ? = 0.01. We have arbitrarily 
chosen the point k = f for the graphical presentation. 


where x,y G M. and 0 < Kg < 1. If we now define the 
iteration of a function for some a: G K as Xq = x, Xi = 
g{x), ■ ■ ■ ,Xn = g(xn-i), we immediately get that 


|x„ - x„+i| < Ar”|xo - xij. (B.2) 

Furthermore, this tells us that for any positive integers p 
and n 


^n+p 


p-1 

m—O 


^n+m ^n+m+l| 


p -1 

<Y^K-+^\x3-X3\ 

m=0 


(1 - ATp 

I- Kg 


jxo - Xij, 


(B.3) 


where we in the first step applied the triangle inequality. 
Since this is valid for any p, we can take the limit p —>■ 00 , 
resulting in 

A" 

|a;„-x*|<-—^jxo-xij. (B.4) 

1- Kg 

Here x* is the so called fixed point defined by x* = (?(x*). 
For contractions, this fixed point is known to exist and 
be unique according to the Banach fixed point theorem. 


If we now introduce the dimensionless parameters e = 
A/A and y = a/^o, we can rewrite the equation for the 
energy of a fixed point in reciprocal space as 

£ = /(x\/l - = ^(e)- 


\9{x)-g{y)\ < Kg\x-y\ 


(B.5) 
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Figure 8. Topological phase diagram obtained by calculating 
the winding number using the ^ ~ approach. The white 
line plots the border between regions where the topological in¬ 
variant is of different parity, as calculated using Eq. (14l; this 
equation is valid for arbitrary coherence lengths. Parameters 
used are an = 0.01, = 50a. 


Evidently, the energy satisfying this equation is the fixed 
point of h. If h is then also a contraction, we can ap¬ 
ply the machinery developed in the previous paragraph. 
Through the mean-value theorem, we know that for any 
a; < z G [ 0 , 1 ) there exists a ?/ G [cc, zj such that 


\h{x)-h{z)\ = \x-z\x 


|y| df{x) 




c=X\/^-l 


(B. 6 ) 


Comparing to Eq. (B.l), then, if 


Kh=X max 


\y\ df{x) 


pi \ Vi — 




dx 


c=xx'i--y 


< 1 , 


(B.7) 

h becomes a contraction. Usually Kh is of the order of 
10 “^ — 10 “^ regardless of the chosen k point, but due to 
the complicated form of / this has to be checked for each 
case separately. Nevertheless, we can finally conclude 


that 


ko-e| < -£i|, (B. 8 ) 

where Eq = h{0), Si = h{so), and e = h{e) is the actual 
normalised energy. Note that the approximation 
corresponds to e —>■ £o- An example of the effect of this 
approximation is seen in Fig. a), where the first ap¬ 
proximation eo{k) and the fixed-point solution s{k) are 
plotted in the same figure. The iterations converge to¬ 
wards the exact solution very rapidly, but, as seen in the 
figure, the first approximation using is already rea¬ 
sonably close. A comparison of the differences between 


iterations for a selected k point is seen in subfigure b); 
the initial error is minor, and each iteration reduces it by 
several orders of magnitude. As expected the errors are 
lowest at high coherence lengths; however, even at low 
coherence lengths, error for the first iteration remains 
lower than 10“^A. 

Another way to test the accuracy of our approxima¬ 
tion, is to compare the topological phase diagram calcu¬ 
lated using Eq. (131 to the phase diagram of the Z 2 invari¬ 
ant found using Eq. (14|. Note that the former method 
relies on the approximation, while the latter is inher¬ 
ently equally valid for any coherence length. As even the 
quantitative differences between iterations are generally 
negligible, the qualitative effects on the topological prop¬ 
erties of the system are minor. In Fig. we see these two 
methods in superposition. The regions in which the two 
methods disagree coincide with areas where the gap size 
is extremely low (see Fig. [^a)). The apparent disagree¬ 
ment between the two phase diagrams is hence due to 
unreliability in the numerical calculation of the winding 
number, rather than the ^0 approximation. Physically 
this is of little relevance, as the gap size in that are is 
too low for observation of topological effects of the type 
studied here. 


Appendix C: The Winding Number formula 


The topological invariant for the one-dimensional FM 
Shiba chain can be calculated using Eq. (12|. This is 


straightforward for the linearised two-band model, but 
the four-band NLEVP is not formulated in terms of a 
Hamiltonian operator. As mentioned in the main text, 
this can be circumvented by defining a topologically 
equivalent flat-band version of the effective Harniltonian, 
H = P+ — P_. Inserting this expression in Eq. (12) gives 




niz ja 


dfc tr 


— •nja 


C(P+-P_H(P+-P_) 


(c.i) 

We can progress by considering the chiral symmetry op¬ 
erator C. The effect of C on the energy eigenstates is 
known to be C\E±) = [Pip). It is straightforward to see 
that this also satisfies CH = —HC. Together with the 
properties of the trace, we obtain 




f'TT ja 


dk tr 


— tt /a 


CP+dkP+- P+CdkP+ . (C. 2 ) 


Integration by parts finally yields Eq. (13) from the main 
text. 
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